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This paper investigates bright quantum-matter- wave soli tons beyond the Gross Pitaevskii equation (GPE). As 
proposals for interferometry and creating nonlocal quantum superpositions have been formed, it has become 
necessary to investigate effects not present in mean field models. We investigate the effect of harmonic confine- 
ment on the internal degrees of freedom, as the ratio of zero point harmonic oscillator length to classical soliton 
length, for different numbers of atoms. We derive a first-order energy correction for the addition of a harmonic 
potential to the many-body wavefunction and use this to create a variational technique based on energy mini- 
mization of this wavefunction for an arbitrary number of atoms, and include numerics based on diagonalization 
of the Hamiltonian in a basis of harmonic oscillator Fock states. Finally we compare agreement between a 
Hartree product ground state and the Bethe ansatz solution with a Gaussian envelope localizing the center of 
mass and show a region of good agreement. 
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I. INTRODUCTION 

Since the experimental realization of Bose-Einstein con- 
densation (BEC) with dilute atomic gases [I, 2], much 
progress has been made in the degree of control possible 
in terms of external potentials and control over interactions 
via magnetic and optical Feshbach resonances. Recently, it 
has become possible to create condensates of atomic species 
with scattering lengths that can be tuned to be negative [ - 
5], for example hyperfine levels of **^Rb, ^Li and '^^Cs [ , ]. 
Chromium is also shown to have negatively tunable scattering 
lengths, but also has significant long-range dipole forces [8]. 

Such attractive condensates have the remarkable property 
of 'self trapping', i.e. being localized (at least in terms of pair 
coiTelations) on a length scale shorter than would be expected 
for a non-interacting condensate. In the limit in which there is 
one dimension where the non-interacting ground state would 
be infinitely wide, the Gross Pitaevskii equation (GPE) pre- 
dicts the ground state would still be localized to a finite size. 
Therefore with the GPE, these condensates behave as solitary 
matter waves, or in the quasi ID case, as classical solitons. 
The transitional region between ID and 3D has been inves- 
tigated using variational methods [9]. Experiments measured 
systems with a lifetime of several seconds, both for the case of 
a single wavepacket [ ] (a ground state) and multiple smaller 
wavepackets [4, 10], refeiTed to as soliton trains. In addition 
to experimental results, theoretical results relating to the inter- 
action of such systems with potential barriers predicted effects 
such as enhanced reflection and transmission [11, 12]. 

In this situation, the GPE predicts an infinite number of 
conserved quantities within the system and thus complete in- 
tegrability [ ]. As a result the inverse scatting transform 
can be used to obtain solutions [_ ] that are a combination 
of bright solitons, which do not change shape as the quantum 
pressure is exactly balanced by the nonlinear interaction, and 
radiation, which does. In fact any initial condition for the GPE 
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equation can be broken up into these components [ ' ']. In the 
case of systems of a multiple soliton system, individual soli- 
tons can collide with other solitons without a transfer of en- 
ergy between them, resulting in only an asymptotic position 
and phase shift. Such localized matter waves (which are typ- 
ically of the order of a few micrometers in width) could also 
theoretically be split coherently into multiple parts [ ' ] and 
prove useful for interferometry [16, 17], studies of quantum 
reflection [12, 18] or probing surface potentials [19]. 

Under certain circumstances such solitons behave like clas- 
sical particles with finite range interactions [ ] even in the 
presence of harmonic confinement. However, in reality, such 
objects should behave as quantum particles (i.e. with no sub- 
structure, but with the location determined by a wave-function 
obeying quantum mechanical laws rather than a specific po- 
sition). Despite the GPEs success in describing many phe- 
nomena in BEC, even for very small numbers of atoms [21], 
this quantum mechanical center-of-mass behavior is totally 
lost under the approximation of a product state wave-function. 
Therefore we consider a full many-body quantum description, 
making use of the usual pseudo-potential approximation. 

The dynamics of the center of mass of an interacting gas 
in a harmonic potential are independent of the interactions, 
giving rise to the so called "Kohn mode" [22]. More gener- 
ally, any potential which looks locally harmonic on the length 
scales dictated by the internal degrees of freedom (in this case, 
the classical soliton length) can be considered to only weakly 
couple the center-of-mass to other degrees of freedom. As a 
result, this behavior needs to be considered separately. In this 
weak coupling approximation the center-of-mass behaves like 
a non-interacting particle of mass NM, which can be localized 
on scales far wider than a classical soliton length or even not at 
all. This delocalisation is not present in the mean field model, 
however the Kohn mode is still present for translational mo- 
tion of the center of mass. 

Exact results exist for Bose gases in free space with peri- 
odic boundary conditions for both repulsive [ ] and attrac- 
tive [24, 25] interactions; these have the advantage of ex- 
plicitly separating the center-of-mass component of the wave- 



function and being accurate even for very small numbers of 
atoms and fragmented states, for which the GPE is not. It 
has however been shown that in certain situations, in the limit 
A^ — > oo, ^ — > with Ng - constant, the GPE functional is 
an exact description of the system [ ]. Modern numerical 
approaches are available, such as the multi-configurational- 
time-dependent Hartree method, which have been used to 
study bright matter wave solitons [ ]. The separation prop- 
erties of the center-of-mass in free space and harmonic traps 
make available the possibilities of creating non local superpo- 
sitions [27, 28] (cf. [ ]) with sufficient numbers of atoms to 
make them detectable, in ways not predicted in the classical 
field description of the GPE. 

In the following we study the Lieb-Liniger(-McGuire) 
gas [23, 24], a ID system of identical Bosons with attractive 
contact interactions, with the addition of a harmonic trapping 
potential. We have present a series of analytic and numeri- 
cal techniques to study many-body effects, making use of the 
separability of the many-body wave-function into a center-of- 
mass component and a relative component. This paper fo- 
cuses on the ground state of the system and energy corrections 
to the relative components of the wave-function as trapping is 
increased, along with estimates of the overlap with the free- 
space relative ground state. 

We consider a unit system with the Hartree soliton length 
set to unity and -\fy a dimensionless ratio between this length 
and the harmonic oscillator length in the axial direction, such 
that y — > recovers the free case solved by the Bethe Ansatz. 
We use a numerical method based on exact diagonalization of 
the Hamiltonian over a basis set of Hermite functions, trun- 
cated up to a maximum energy, in order to determine the 
many body ground state energies, combined with a variational 
method for low y. 

This paper is organized as follows: Section II introduces 
the exact results in one dimensional infinite systems (using the 
Lieb-Liniger model [ ]) and the unit rescaling used to keep 
the mean field soliton length constant throughout the paper 
Also included is the separability of the many-body Hamilto- 
nian and the existence of the Kohn mode as well as the exact 
eigenstates for two interacting bosons in a harmonic potential. 
Section III derives a perturbative energy correction to the rel- 
ative ground state energy from the introduction of a harmonic 
trapping potential, along with a variational procedure to esti- 
mate the ground state in the limit of weak trapping. Section IV 
introduces the numerical method used to perform calculations 
in the many -body system for varying ID harmonic trapping 
potential, using a basis set of harmonic oscillator eigenstates, 
which are projected to a center-of-mass excitation basis. Sec- 
tion V numerically investigates changes to relative component 
(i.e. having excluded the center of mass) of the ground state 
as the trapping potential is increased. These calculations are 
performed for different numbers of atoms and compared with 
predictions based on the GPE. Section V C examines quantita- 
tively the overlap between the mean field approximation of the 
ground state and the many-body solutions, along with finding 
a regime of agreement where the difference between the two 
models is small in every respect. Section VI summarizes and 
comments on the results. 



II. PRELIMINARIES 
A. System overview and rescaling 

1. The many-body Hamiltonian in first and second quantization 

We consider a system of identical bosonic atoms within 
a cylindrically symmetric prolate (the radial frequency o),- is 
greater than the axial frequency Wv) harmonic trapping poten- 
tial V(x,y,z) - M[aj]x^+ajj(y^+z^)]/2, where Mis the atomic 
mass. We further consider the system to be sufficiently low- 
temperature for the atomic interactions to be pure i-wave and 
contact-like, and assume that we are in an appropriate param- 
eter regime so that the radial modes can be considered "frozen 
out" for low -energy states, taking the Gaussian form of radial 
harmonic oscillator ground states [30-32]. Finally, we assume 
the interactions to be attractive. 

Integrating out the radial degrees of freedom, we obtain the 
ID Hamiltonian, in second-quantized form, and hence for an 
arbitrary number of atoms, as follows: 



-/ 



dx^'\x) 



'iMdx^ 



MljIx^ 



^1 



^-ix) 



£/jc>I'f(jc)»I'f(x)>P(x)>P(x), (I) 



where gi£, - 2TicOr\as\, and a^ is the (assumed negative) .s-wave 
scattering length [ ]. The coordinate-space representation 
for the corresponding first-quantized form of this Hamiltonian 
for A^ atoms is then 



m.?) = £ 



k=\ 



'2MaS| 



MojWk 



N k-1 
■ glD ^ ^ S(Xk - Xj), 
k=2 j=\ 



(2) 

where we use x as a shorthand for {jci, X2, x^, . . . , x^}, the A^ 
individual particle coordinates. In the absence of any trapping 
potential (i.e., w^ - 0), the system is formally integrable [ ], 
and the first exact results for the case of attractive interactions 
were obtained in 1964 by McGuire [24]. 



2. Hartree factorization: The Gross-Pitaevskii equation 

In this approximation one assumes the many-body wave- 
function ifr(x) to be factorizable into product form, such 
that each individual atom is described by the same single- 
particle wavefunction (p(xk). Hence, the Hartree wavefunc- 
tion tffH(x) - Y\^=i (pixk)- Minimizing the energy E - 
J dx((r'^(x)H(x)if/ii(x) of such a stationary state with respect 
to variations in 0(xi) leads to [34] 



jt/0(x) = 



n^ 



2v2 



'iMdx^ "^ 2 



;\D 



(N-imx)\' 



4>(x\ 



(3) 



with (p{x) normalized to unity, and // a Lagrange multiplier, 
given by 



/ 



dx(f>*(x) 



'iMdx^ 



Mtolx^ 

2 

•^id(A^- 



P(x) 
l)jdx\c/>(x) 



(4) 



Equation (3) is the one-dimensional time-independent Gross- 
Pitaevskii equation (GPE) [ ], which formally tends to an 
exact description as A^ — > oo while gmN is held constant 
[35, 36]. In this limit gioiN - 1) « gioN, and it is more 
typical for the coefficient of the nonlinearity in Eq. (3) to be 
set proportional to A^. As we will also consider small particle 
numbers, in what follows we choose to retain the proportion- 
ahty to (A^- 1). 



3. Resettling to dimensionless form 

It is convenient to rescale our description of the system in 
terms of an effective K - M - gmiN - 1) = 1 unit system, 
referred to as "soliton units" [20, 30]. Space, time and energy 
scales are then given in units of fp-/Mgij:,{N- 1) (the classical 
soliton length [ ]), Tr' IMg\^{N - 1)^, and Mg\^(]^ - iflh^, 
respectively. 

We work within this system of units from this point on- 
wards. Equation (1) then simplifies to 



H 



/***« 
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2dx^ 2 

— - \ dx^\xyvHxy^{x)mx), (5) 



2(A^ 



the first-quantized form of the Hamiltonian [Eq. (2)] trans- 
forms to 



Hi^^Yu 



k=\ 



\_(f_ fxl 

'idx^^ 2 



N k-l 



-]^ZZ^(^^-^^(6) 



k=2 j=\ 



and the GPE [Eq. (3)] becomes 



'2~dx^ 



2 2 

y X 



\4,{xt 



4,{x). 



(7) 



We have introduced the dimensionless parameter y, which 
is the square of the ratio of the classical soliton length to the 
harmonic length yjtilMoj^ [30], i.e.. 



r = 



Mg\^{N - 1)2 



(8) 



Within our chosen system of units y appears in the rescaled 
Hamiltonian and GPE as a dimensionless effective trap fre- 
quency. This also reveals y to be the only free parameter in 
the GPE, which as a description of the system is effectively a 
classical field limit, and the particle number A^ appears as an 
additional free parameter in the fully quantal Hamiltonian. 



B. Known exact results 

1. Exttct results in free space 

In the case where there is no axial trapping potential, i.e., 
y - Q, the many-body eigenstates [ , ] and GPE station- 
ary states [ 1 4, 30, 34] are known. The stationary solutions to 
the GPE [Eq. (7)] minimizing the energy are classical bright 
solitons [_ , _ ] 

0(x)=isech(^), (9) 

where the value of xq is arbitrary (as is that of an irrelevant 
global phase). The Hartree approximation to the many-body 
wavefunction is thus 



H(f).i,nsech(^), (10) 



«A: 



and is localized around xq, i.e., i^h(-^ ^ as \xk - xo\ — > 
oo. Exact solutions also exist in box and periodic boundary 
conditions, given by Jacobi elliptic functions [ ]. 

The exact ground state wavefunction for Eq. (6) with y - 
(see Appendices A and B 2) is an A^ -particle bound state, 
proportional to [24, 37] 



</'g(-^ 



{N-\)\ 



\N-\y 



nexp 



N k-\ , ,\ 

Vn Vn \Xk — Xj\ 



(11) 



As the Hamiltonian [Eq. (6)] has no external potential, all its 
eigenfunctions are independent of the center-of-mass coordi- 
nate 



xc 



1 ^ 



(12) 



(see section 11 B 2) and there exists a continuum of 
moving A^ -particle bound state eigenfunctions iJ/{P,x) = 
e'''''^if/G(x)/ y/lK. The normalization convention is then such 
that ] dxip*{P, x)i//(P', x) - 6{P - P'), and the ground state is 

written as i^(0, x) = iffcix)/ y/ln. The exact ground state is 
thus completely delocalized in the center-of-mass, and is only 
localized in the sense that tf^cix) — > as \xii - Xj\ — » oo, for 
any k, j. 

The localized Hartree solution (Ah(^ violates the transla- 
tional symmetry requirement imposed by the absence of exter- 
nal potentials in the Hamiltonian due to the tacit assumption 
that the minimizing wavefunction should vanish as xj — ^ ±oo 
[ ]. Note, however, that the particle densities about a spec- 
ified value R of the center-of-mass coordinate [given by the 
expectation value of 6{R-xc) 2^i d(x-Xk)] corresponding to 
tffuix) [Eq. (10)] and ifrdx) [Eq. (11)] agree to order 1/A^ [.■ .], 
and hence are identical in the limit N ^ oo. The energies En 
and Ec corresponding to the wavefunctions i/'h(^ and if/cix) 
are given by 

N 
Eh^-wJ' (13) 



24 

A^(A^ + 1) _ N 

rLn — = £- H ~ ■ 

24(A^-l) 12(N-1) 



(14) 



As one would expect, the exact eigenenergy Ec is less than 
En, and the difference in energies per particle (En - Eg)/N - 
1/12(N - 1) vanishes as A^ — > oo. 



2. Separation of the center-of-mass coordinate Xq 

In the case of any external potential being either harmonic 
or nonexistent, the center-of-mass dynamics separate and are 
independent of any two-body interactions. Consequently 
the center-of-mass eigenstates are simple harmonic oscillator 
eigenstates or plane waves, respectively, in the former case 
this is referred to as the Kohn mode. This may be readily 
seen by expressing the first-quantized form of the Hamiltonian 
Eq. (6) in terms of the Jacobi coordinates, i.e., xq [Eq. (12)] 
together with 



1 *"' 

7=1 



(15) 



for k G {2, 3, 4, . . . , A^). The Hamiltonian can then be phrased 
as i/ = He + ^R, where 



1 52 A^r'^c 



(16) 



N 



HR(i)^j] 



k=2 



k d^ (k- Dy^e, 



2(k-l)df, 



2k 



N 



1 N ( k-l ^ ^ 



(17) 



A' 



k=2 
. N k-\ ( k-\ ^ . . \ 

k=2 j=2 V {=j+l ■' 



^ is a shorthand for {^2, ^3,^4, ■■ ■, ^n), and we have used the 
identity x, - xj = ^, + Y^tU ^(je - [(j - ^)IMj (with b>a 
and ^i = xq). In cases where the upper limit of a sum is less 
than its lower limit, the sum is taken - 0. 

Hence, the normalized ground state of He is exactly 



iAc(.xc) = ( ^^ 1 exp I - 



(?) 



Nyxl'^ 



(18) 



with eigenenergy = 7/2. 



3. Two interacting bosons in a harmonic potential 

The case of two identical bosons in a harmonic potential 
with contact (5-function) interactions is also exactly solvable 
[40, 41]. In this case the eigenfunctions of Hr(^2), defined 
through H^{^2)M^2) = £'r,«0«(^2), are given by 



P„(^2) = N„U(-y„, l/2,r^2/2)e-rf?/4^ 



(19) 



implicit solutions of 



r(l/2-y„) 



1 



r(-v„) 2^27' 

and set the eigenvalues of //r(^2) through 
Er,„ = 2v„ + 2 p- 



(20) 



(21) 



Attractive interactions must reduce Es,q from the noninteract- 
ing case, so that Er^ < yjl => vq < 0. As outlined in ap- 
pendix C, it then follows that in the limit y — > (interaction 
dominated regime) Erq — > -1/4 H- 0{y^). This is in agree- 
ment with the total ground state energy Eq for the case of 
two attractively interacting bosons in free space [Eq. (14)], 
as one would expect due to the center-of-mass energy of the 
free space ground state being = 0. In the opposite limit of 
y ' — » (trap dominated regime) harmonic oscillator eigen- 



values and eigenfunctions must result, i.e., E„ 

e2, 



{ivi + i/2)r 



and U{-v„,\l2,yfJ2) -^ //2„( Vr^2)/22"/ V2, ^j^g^g ^^ 
H2n are even Hermite polynomials.' 



ni. PERTURBATIVE AND VARIATIONAL METHODS 
A. Interaction dominated limit in a harmonic potential 

In the case where 7 «; 1, we may consider the effect of 
the trap to be dominated by the effect of the interactions, and 
therefore negligible in //r. As there are no interactions present 
in He, the effect of the trap is in this case always significant, 
even in the interaction dominated regime. 

We may therefore consider a limiting case Hamiltonian Hq, 
composed of //r [Eq. (17)] with y - Q, plus He [Eq. (16)]. 
Written in terms of conventional single particle coordinates. 



Hq{^ 



1 viL 



r 



V 



l^x, 



A=l 



1 



A'- 



N k-l 

tZZ 

k=2 j=l 



5{xk-Xj), 

(22) 

and the correctly normalized ground state i/^o can be put to- 
gether from Eq. (11) multiplied by Eq. (18), i.e., t/ry = (/'c'/'G; 
with the sum of the corresponding eigenvalues determining 
the overall energy £0. Hence, in terms of single particle coor- 
dinates 

( r « n2\ 



\I/q{x) 



(?) 



1/4 



exp 



'2N 



Z 

k=\ 



Xk 



^g(A (23) 



and, from Eq. (14) plus y/2 (the harmonic oscillator zero 
point energy), 

N(N + 1) 



r 



(24) 



24(N - 1) 2 

This interaction dominated limit does not correspond to any 
physical system but is a useful starting point for perturbation 
theory. 



where U(a,b,z) is the Tricomi confluent hypergeometric 
function [~'^], and N„ is a normalization constant. The v„ are 



Due to Bose symmetry <fy„(^2) = <Pii(-^2), ie-, eigenfunctions must be even. 



B. Perturbation results 

To proceed from the approximated Hamiltonian (23), we in- 
clude the effect of the harmonic trap on the relative degrees of 
freedom via Rayleigh-Schrodinger perturbation theory. The 
full Hamiltonian (2) can be written as H{x) - Ho(x) + AH(x), 
with 



AH(x)^ 



N ] ( ^ \ 

k=l \k=l J 



(25) 



As AH{x) oc y2, we expect perturbation theory to yield partic- 
ularly good results in the limit of small y. For the first-order 
energy correction to the ground state, 



E'" = mAH\^o) 



I 



dx iJ/q(x)* AH(x)iJ/q{x) 



(26) 



which also serves as a definition of the bra-ket notation, we 
find (Appendix B 4): 



£(!> = y 



2(A^_-_1) 



2 W-l 



§^^' 



(27) 



The sum in Eq. (27) is simply the second Harmonic number, 
for which the asymptotic behavior in the N » I limit is given 
by [4.J] 



4-^k^ 6 N-\ ^ ' 



Thus, asymptotically the energy correction goes as 



(28) 



£(!> ~ / 



>-T-i-^M 



(29) 



For large A^, this coincides with the result obtained using the 
free space Hartree solution, given in Eq. (10), as an approxi- 
mation for the ground state (Appendix D) 



7(1) 



(A^-1) 



f 



2 ,,2 ^2 



dx 



sech(x/2) y X- 



= /(A^-l)- 



(30) 



These results are displayed in Fig. 1; for small A^ there is a 
large difference between the result predicted by the Hartree 
product state [Eq. (27)] and the result predicted by the exact 
many-body ground state [Eq. (30)] of the approximate Hamil- 
tonian (22). There is also a weak number dependence from the 
Harmonic series in Eq. (27). However as A^ » 1 both meth- 
ods give the same energy correction per atom, jp'y^ 16. For 
A^ = 1000 the relative difference (£[]' - £'(")/£'(i> ^ 0.0016 is 
already small. 




-Many-body first order energy correction 
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FIG. 1. (Color online): First-order energy correction per atom for a 
many-body soliton with the center-of-mass in the lowest eigenstate 
of a harmonic oscillator given in Eq. (23). The exact solution given 
in Eq. (27) and its expansion up to next-to-leading order, given in 
Eq. (29), begin to agree for A' ^ 10 with the latter always underesti- 
mating the true value. Both curves approach the approximate result 
predicted by the Hartree approximation [Eq. (30)]. The relative dif- 
ference between different predictions lies below 1% for A' ^ 165. 



C. Variational minimization 

In order to improve the value for the ground state energy 
beyond the first-order-perturbation-theory result (27), we use 
the (normalized) variational ansatz 



r.i(^ = Mxc)NA- 



(N-l)/2 



exp 



V \<k<j<N 



\Xk 



2[N-l] 



with y > and the constant N the same as Eq. (11)^, 



Ar = 



N-\ 



'{N-D 



(31) 



(32) 



which is calculated in appendix B 2. Since the center-of-mass 
wavefunction is unchanged, we will only have a correction to 
the relative energies, these are calculated in Appendices B 3 
and B 4. The total energy for this wavefunction is 



<iAti;.|//|^tl>) = [lA -A^)Eg + ^ + \ 



(33) 



for the expectation values of each section of the relative 
Hamiltonian, with Eq being the (negative) ground state en- 
ergy of the free soliton in soliton units given in Eq. (14) and 
£''^'* the first-order correction given by Eq. (27). In order to 



^ With just N as a prefactor the equation would be normalized with respect 
to Ax, scaling this out along with y — > yjA^ to keep the center of mass the 



same, gives the extra factor of /t*'^ '"^ 



calculate the energy minimum, the derivative of Eq. (33) with 
respect to A has to be zero: 



(2 - 2A)Ec - -^ = 



(34) 



which (for /I t^ 0) is equivalent to a 4th order polynomial in A 

A'^-A^ -K^O, (35) 

where the constant k is defined as the ratio of the ground state 
energy and first-order correction 

224(N-l) 



7 



(N + 1)N^ 






(36) 



For fixed N, k ccy^, the value of this prefactor is an increasing 
function of A^ with a minimum of k = 2y^ at N - 2 with an 
asymptotic limit of [cf. Eq. (29)]: 



lim K - 4-t: y 



(37) 



Thus, K is small for y <K 1. 

Equation (35) has four roots, only one of which is real and 
positive, which is the root of interest. The exact analytic solu- 
tion is given in Appendix E. 

If K <^ 1 , this solution is approximately 



10 



1 +K 



(38) 



Eq. (40) still is A^-dependent due to the A^-dependence of Aq 
[Eq. (38), cf. Appendix E]. 

Figure 2 (a) shows the overlap (40) as a function of y for 
various particle numbers. For /lo, the exact value given in Ap- 
pendix E was used. As expected, the A^-dependence is quite 
strong. Figure 2 (b) shows the A^th root of the overlap (40), 
i.e. the effective single-particle overlap. The effective single- 
particle overlap is larger than 0.99 f or y < 0.15 for all A^, in- 
dicating that i/'o(-*t) from (23) is still a good description in this 
parameter regime and the trap has had little effect on the inter- 
nal degrees of freedom. The limit A^ ^ oo is nearly reached 
for particle numbers as low asN = 100 [note that in panel (a), 
the limit A^ — > oo would lie on the coordinate axes]. 



which leads to the minimum in the energy of 

„2 (£(1))2 



(1) 



2 Eg 



+ 0([E'^^r/\Ecn. (39) 



As the variational Ansatz (31) does not affect the center- 
of-mass part of the wavefunction, calculating the overlap be- 
tween this variational Ansatz and the state Eq. (23) (i.e. the 
A - I state) is an interesting physical quantity: its modulus 
squared is the fraction of the relative wavefunction which is 
projected to the relative ground state if the trapping potential 
was turned off quasi-instantaneously (cf. [-n]). The overlap is 
given by (see Appendix B 3): 



(•AlUVo) = 



Al''-+A, 



1/2 



N-1 



(40) 



Using the approximation (38), the overlap (40) approxi- 
mately is [1 + K- /8 +0{K^)y^^'^^ and we thus expect the overlap 
to vanish in the limit A^ — > oo for k > Q. Rather than inves- 
tigating the total wavefunction overlap (40), the A^th root of 
Eq. (40), an effective single particle overlap, is a more suitable 
value in the limit A^ » 1 as it tends to a constant as A^ — > oo 
and is related to comparing two GPE orbitals. Note that for 
two Hartree-product wavefunctions, the effective single par- 
ticle overlap would be independent of A^, but the A^th root of 




FIG. 2. (Color online): (a) Total wavefunction overlap, given by 
Eq. (40), (b) effective single-particle overlap, given by the A'th root of 
Eq. (40), of the variationally obtained solutions for different rescaled 
trap frequencies y with the free space ground state solution (y = 0) 
with a Gaussian envelope for the center of mass. Effective single- 
particle overlap is treated as the A'th root of the total overlap as for 
two different product states this is independent of number and equal 
to the overlap between the single particle wavefunctions. Bottom to 
top the solid lines on both graphs correspond lo N = 100, 10, 6, 3, 2, 
the dashed line corresponds to the A' — > oo limit of the variational 
many-body solution [using k from Eq. (37)] and is very close to the 
N = 100 line. 



work well for y > 0.16. 



B. Computation procedure 

We expand the field operator over the set of Hermite func- 
tions 



n(Wx) = 



w 



;t!2%i/2 



Hk(Wx)exp[-W^x^/2) , (41) 



where i/^ are the Hermite polynomials, giving ^'(x) - 
Yjk fk{x)ak- The Hamiltonian (5) in this basis can be split into 
three separate parts: 

W^ v-i r t 



- yj(k + l)(k + 2){al^^ak + fl|fli+2)] 
the kinetic Hamiltonian, 



t<^k 



k 

+ yl(k+ \)(k + 2)(a[^^ak + fl^flt+i)] 



(42) 



(43) 



IV. COMPUTATIONAL METHODS INCLUDING A 
HARMONIC POTENTIAL 

A. Overview 

While the focus of the previous section lies on the case 
of small y, the numerical methods introduced in this section 



the potential Hamiltonian, and 



W 



k£mn 



(44) 



the interaction Hamiltonian. The factor of fkimn is the integral 
of four Hermite functions (with W set to unity) over all space, 
i.e. 



-! 



fktmn = dx ipkix)(pe(x)ip„(x)ip„(x) 



(45) 



Here, the functions are also real, so there is no need to take the complex conjugates. This can be calculated exactly in terms 
of gamma functions r(x) and a standard hypergeometric function 3F2 evaluated at unity [ ..] 



/* 



1 



ktmn 



-T{[k + (-m + n + l]/2)r([/t -( + m-n+ \]l2)T([-k + l + m-n+ l]/2) 



^|2^J:^ ^ k\{\n\(m - n)\ 

X 3F2([-n, (m-n + k-e+ l)/2, (m - n - k + £ + l)/2]; [1 + m - n,(m - n - k - { + l)/2], 1) . (46) 



I 
Without interactions, the ideal gas Hamiltonian is given by be expressed as 

Hideal = i/K + Hp. (47) 



For W = -y/y, the basis states (23) are eigenstates of the non- 
interacting Hamiltonian. The total Hamiltonian can therefore 



^ = rZ(*4)' 



fl^flj. - 



Vr 



2(A^-l) 



/ ^ fkemn^jfifim^n 7 (48) 



k^mn 



and we refer the ground state of this as 1^^2(7)) and the ground 
state energy as 



{iff.iyWmy)) = Egiy) . 



(49) 



C. Truncation and projection to center-of-mass excitation 
basis 

In order to do computations we must only use a finite ba- 
sis set, which will introduce the inaccuracy. This is discussed 
in Appendix F. Essentially all possible states for which the 
eigenenergy related to the Hamiltonian (47) lies below an en- 
ergy cut-off Ecut are included, we refer to this set as the 'trun- 
cated basis'. 

The brute force procedure would now be to calculate all 
the matrix elements of the interaction Hamiltonian (44) using 
this truncated basis and add the matrix of energies of the ki- 
netic and potential Hamiltonians and diagonalize this to get 
the eigenstates and energies. However we can reduce the size 
of the truncated basis set used in this computation by recall- 
ing from Eq. (16) that the center-of-mass Hamiltonian com- 
mutes with the relative Hamiltonian and thus they have sepa- 
rate eigenstates. 

Inspired by the ladder operator treatment of a single particle 
in a harmonic oscillator, we define 



A^ixc) - \\ \Nyxc + - — 

V 2Ny \ ^ dxc 



(50) 



Noting we can express the center of mass Hamiltonian as 
Hc(xc) = 7(A^ixc)A^{xc) + 1/2), as is the case of the single 
particle ladder functions. Moving to second quantization, we 
can equivalent operators in terms of creation and annihilation 
operators in our basis (for W = -\fy) via [ ] 



= Xi ^^+^«l«t+i 



(51) 



where A+ - (A")"'". These satisfy [//o,A*] - +7A* and thus 
energy levels spaced in units of y; they also commute with 
the interaction Hamiltonian [Hj,A'^] = as required. Note 
that for A^ - 1, A+ acting on the ground state can be used 
to construct all the eigenstates of the system. We then can 
express the center-of-mass Hamiltonian as 



He 



7 

N' 



A^A- 



7 



(52) 



In order to use this property to reduce the basis, we must first 
transform our truncated basis (E < E^ui) to basis states which 
are eigenstates of the center-of-mass Hamiltonian, this proce- 
dure is detailed in Appendix F 3 along with the reduction in 
basis size it achieves. 

We keep only eigenstates of He with eigenvalue y/2 
(center-of-mass ground state) and project the Hamiltonian 
(48) from the truncated basis to this new and reduced ba- 
sis. For high A^ > Ecut/7 - N/2, the reduction asymptotes 
to 7t/ -y/6(£'cut/r - N/2) for E^^Jy - N/2 » 1. For details see 
Appendix F 4. 



D. Using dlfTerent-wldth Hermlte functions 

Using functions with a.W - ^y, such that they are eigen- 
states of Hjdeab is not desirable in the y — > limit because 
the basis will consist of states much wider than the wavefunc- 
tion we are using them to construct. For an infinite basis, the 
ground state should be independent of the basis used to de- 
scribe the system (in our case, it should be independent of the 
value of W). For numerical calculations, the basis will be fi- 
nite and thus some choices of W are better than others. In 
Sec. V, we will calculate the ground state for y = in order 
to determine the optimal value for W to be used in the calcu- 
lations. 

For arbitrary W, the Hamiltonian now reads; 



«=z 



W^ + 7^W-^ 



k + - \aj^ak 



7^W-^ - W- 



^J{k + l){k + 2){a[ cik + Sifli+a) 



W 



2{N-\) 



rr^Z-^* 



Mmn^j^^l ^m^n » 



(53) 



klmn 



which includes extra mixing terms in the ideal gas Hamilto- 
nian (47). This causes a fairly significant issue in that it is no 
longer possible to exactly separate center-of-mass eigenstates 
in this basis, meaning the full basis would need to be used in 
order to achieve the center-of-mass ground state. Using this 
method with just a truncated Hilbert space and no projection 
to the sub space with zero center-of-mass excitation would 
make achieving convergence painfully slow. The solution to 
this is therefore to reduce the basis in the same way as be- 
fore, but accept that the center-of-mass wavefunction we end 
up with is given by 



fcixc) = 



W 



exp 



-A^- 



W 



l.^^ 



(54) 



which is not an eigenstate and has energy 
Eq - {W^ + 7^W^^)/A rather than the true y/2, thus we 
know the true ground state is the wavefunction we obtained, 

multiplied by J ^Jy/W &x^{[7 -W'^]Nx^/2\. This approach 
has the huge advantage that, if W is kept constant, the occu- 
pation of the basis states for the ground state should change 
very little as y — > where they will tend to the solutions on 
the infinite line. 



E. Numerical ground states within the GPE approximation 

Within the GPE approximation, we can obtain the ground 
state by solving Eq. (7) as the ground state is the only sta- 
tionary state of the system. The method used here is to again 
expand over a finite basis set of Hermite functions of arbitrary 
width scaling W 



= ^c,Vw.^,(Wx) 



(55) 



*:=0 



then to produce a set of 77 + 1 nonlinear equations in the coef- 
ficient set c (which will be real), by integrating Eq. (7) multi- 
plied by ipkix) over all space, for k - {0,. . .,7]}, giving 

= - jUQ + ^ (k + 1 /2)c, 



y2^-2 - W^ 



( ^|(k^l)(k^Ck^2 + ^k{k - 1)Q_2) 



W 



2{N - 1) 



Z/* 



k^mn^f^m^n^k 



(56) 



{.m,n=0 



There is also an (7/ + 2)th equation, relating to the normaliza- 
tion Yjk kil^ - 1- Denoting the vector with an equation at each 
position as F(c), we wish to solve F = 0. We use Newton's 
method (as in [ ]) to iteratively solve for c, via 



(n)\ /„(«+!) _ „(«) 



Jic'">) (C' 



c^"0 = F(c^"0 



where J is the // + 1 by 77 -H 2 Jacobian matrix associated with 
F. 7] is increased until convergence is achieved. 



V. EFFECTS OF HARMONIC CONFINEMENT 
A. Ground state energy 

Using the methods from the previous two sections, we in- 
vestigate the effect an external potential has on the relative 
component of the ground state |i/'g)(y) [cf. Eq. (49)]. This 
is important to quantify how soliton-like the state is, along 
with what excitations can be expected if the state is released 
quasi-instantaneously from the potential. Such dynamics have 
already been considered using the GPE in [ 4]. 

Figure 3 shows AE/N, the energy difference per atom be- 
tween the numerically calculated ground state energy £'g(y) 
and the ground state energy of the artificial Hamiltonian (23) 
given by Eq. (24), for a range of y and A^ values. It is pro- 
duced by calculating the ground state energy via the three nu- 
merical methods, namely exact diagonalisation in a basis of 
hermite functions with either optimized widths for weak trap- 
ping (shown in table 1), or widths which are eigenstates of the 
non interacting problem, and variational minimization, for a 
range of y and taking the smallest value. This is because, due 
to the variational principle, all of these techniques produce 
only values greater than or equal to the ground state energy, 
hence the lowest is the best estimate. 



B. Universal behaviour 

If we consider instead a rescaling f = -/(N - l)^/N^ and 
AE = AE(N - l)/N (i.e. converting to a unit system in which 
^idA' = 1 as opposed to gioiN - 1) = 1), a more universal 
behaviour is present in AE, with little number dependence as 
shown in fig. 3 (b). To see this analytically, we note that for 
-y <K 1, our variational result of Eq. (39) for the energy is 



A' 


1 


W 


Reduced basis size 


3 


84 


2 


631 


6 


38 


1 


3009 


10 


28 


0.5 


2534 


100 


24 


0.5 


1575 



TABLE I. This table shows the parameters used in calculating the 
graph in Fig. 3, where A' is atom number, 77 is the cut-ofF and W is 
the width taken for the fixed width calculations (coarsely chosen to 
minimize the ground-state energy at 7 = 0). The reduced basis size 
is the number of states (with zero center-of-mass excitation) used in 
the exact diagonalisation of the Hamiltonian. The basis are chosen 
to be a reasonable computational size, however the numerics are less 
reliable for small 7 



applicable. AE is obtained by subtracting the factor of Eg + 



('57") 7 /^^ then converting to our rescaled units we have 



AE J N \r-i 
N ' 



N -\ i-^k^ 

k=l 




-4 2AN^ 




' {N^-\){N-\) 



+o{r) 



(58) 



the A^ dependent factor of order y^ (which is the rescaled first 
order energy correction) is 2 for A^ = 2 and decreases mono- 
tonically to 7r^/6 ~ 1.6 as A^ ^ 00, hence for very small y, 
the N - 2 line is largest, but the difference is very small. The 
order y'^ term has negligible number dependence and so is un- 
likely effect the ordering of these lines within for the range 
of variational models validity. On the other end of the scale, 
as y ^ 00 (the trap dominated system) we can neglect in- 
teractions in Eq. (5), giving a ground state of a product of 
Gaussians of width 1 /y, subtracting the center-of-mass energy 
gives AE — > y(A^ - l)/2N + 0{ yjy) or, in our rescaled units, 
AE — > y/2-i-<9( Vy) and so to leading order, the A^ dependence 
vanishes. 



C. The classical soliton limit 

As shown in Fig. 3, for low y the variational ansatz (31) 
gives the best estimate for the ground state energy of all the 
methods used in this paper For low enough y, this variational 
ansatz in turn is very close to the product of the free many- 
particle solution with a Gaussian center-of-mass wavefunc- 
tion (23) (see Fig. 2). In the limit of small y, the integral 



S 



dx\ ■ ■ • I 

00 %J —Ck 



dxN il/o{x)ifrn(x) 



(59) 



can thus be used to investigate deviations of the Hartree- 
product wavefunction (10) from the true many-particle ground 
state (which is well approximated by iJ/q{x) for y <& 1). As 
was the case for Fig. 2, the effective single particle overlap 
gi/'V will also be considered as we are interested to see how 
well the wavefunction is described by a product state, and 
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FIG. 3. (Color online): Difference between the ground state energy per atom and the energy of a free many body ground state with a Gaussian 
center-of-mass profile, (a) is in terms of EqIN [as defined in Eq. (24)] as a function of rescaled trapping strength y and (b) with a rescaled 
y = y(N - if IN- and A£ = l^E(N - 1)/N . From bottom to top the lines on (a) are N = 2,3, 6, 10, 100 and the dotted top line is the GPE 
prediction (which will agree with the many body results asN ^ oo) outlined in Sec. IV E. The markers indicate a point on the line generated by 
different methods, (red) circles use the variational solution Eq. (33), (green) triangles use the fixed width basis sets to find the lowest eigenvalue 
of the Hamiltonian (53), cf. table I. (Blue) squares are obtained using the basis of eigenstates of the non interacting Hamiltonian (47) to find 
the lowest eigenvalue of Eq. (48). The N = 2 line is plotted using the exact solution [Eq. (21)] detailed in Sec. II B 3 and the mean field line is 
obtained by the method explain in Sec. IV E. The inset shows a zoom of the low y section, demonstrating the initial quadratic dependence on y. 
Figure (b) shows the universal behaviour present using rescaled units, this is the same data as in (a), however the numerical lowest eigenvalues 
are plotted as points to make them visible. 



when comparing two product states with different single par- 
ticle wavefunctions, this quantity is constant with changes to 
number. 

In order to make an educated guess about what range of 
y will give a large overlap, we look at the expectation value 
of the square of the center-of-mass location over the Hartree- 
product wavefunction (Appendix D): 



(^hI^I-Ah) = 



3N 



(60) 



This value is identical to the variance of the center of mass, as 
both the many body and Hartree states are centered about x - 
0. A variance calculation can also be performed for Eq. (23), 
this is particularly simple as the center-of-mass is explicitly 
separate and is given by 



{i//o\xi\il/o} 



2yN 



(61) 



As the Hartree product state is uncorrected, for a large enough 
A^ the distribution associated with the center-of-mass location 
will therefore tend to a Gaussian (with variance of n^/3N) via 
the central limit theorem. We therefore consider an effective 
center-of-mass wavefunction that is the square root of this dis- 



tribution, yielding the overlap integral 



1(7) 



L '"^ ( 



27r4 



exp 



-A^ 



(r + 3/27r^)xl 



(247r>) 



■2-,Al/4 



V2y7r2-H3 
which reaches its maximum / (ymax) = 1 for 

3 



(62) 



Tmax ~ 



27r2 
0.15 



(63) 



As the above analysis focuses on the center-of-mass part of 
the wavefunction, and thus Eq. (62) is likely to overestimate 
the overlap S as defined in Eq. (59), Eq. (62) also predicts 
a y — > behaviour of the form /(y) ~ cy'^"^ with c a con- 
stant. Figure 4 shows a numerical calculation of S for a range 
of y and A^, the integration is performed via Monte Carlo 
methods, i.e. weighted sampling using random variables with 
a sech(.x/2)^/4 distribution (obtained via the ziggurat algo- 
rithm [ ]) until a standard error of < lO"'* was obtained. 

It can be seen from Fig. 4 (a) and (b), that maximum over- 
lap occurs just slightly above y = 0.16 [close to the analytic 
estimate (63)] and improves as A^ increases. Based on our 
previous discussion of effective center-of-mass width, the top 
value should relate to the overlap of the relative degrees of 
freedom, although this is not well defined. Graphs (c) and (d) 
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FIG. 4. (Color online): (a) 2D projection of overlap, S as defined in Eq. (59), between many-body free state with a Gaussian envelope (of 
width oc l/-y) given in Eq. (23) and the mean field soliton solution, given in Eq. (10), for a range of A' and 7. (b) shows Horizontal slices 
through (a), and the dash-dotted line is the analytic estimate, based purely on center-of-mass position uncertainty, given by Eq. (62). (c) shows 
effective single particle overlap S''" and (d) shows the residuals 1 - S''" for given A' values again via slices through (c). The solid lines in 
the lower figure (b) correspond Xo N = 2, 3, 10, 100, 1000 in that order from bottom to top, and this ordering is reversed for figure (d). As 
expected, the effective single particle overlap plot (c) show a rapid convergence to unity as N increases, [(a) and (b)] suggest that most of the y 
dependence in S is due to the efi"ective "center-of-mass width" of the Hartree product solution, since the shape of each overlap curve is similar, 
besides a small offset, to the dash-dotted line. This indicates that a Hartree soliton is a very good approximation to the many body solution if 
the center-of-mass wavefunction is also localized. 



show the Mh root of (a) and (b), effectively overlap at the level 
of single particles, which tends extremely rapidly to unity as 
A^ increases for any y over the range shown. 

A useful point that this high overlap implies is that the 
many-body state Eq. (23), is extremely well approximated 
by the Hartree product state if the center-of-mass envelope 
squared is approximately the statistical distribution that would 
arise from taking the mean of the A^ independent probability 
distributions |0|^ [with (p given in Eq. (9)] associated to sin- 
gle atom positions in the product state. For this reason the 
Hartree product state would be expected to well approximate 
the ground state of the system, even at a many-body level, if 
the center-of-mass envelope is localized to the size of this dis- 
tribution (i.e. y » 3/27r^) by the potential. The converse to 
this is also true, an initial condition that is given by Eq. (10) 
is well approximated by Eq. (23) with y k ijliP'. This could 
be used to estimate center-of-mass position uncertainty of a 
state, initially given by a Hartree product wavefunction, as 
it evolves in time, using known results for the spreading of 



Gaussian wavepackets. 

We also consider how these many -body effects would affect 
experimental observations. From a measurement the atomic 
density, one could use the mean of this signal to determine 
the center-of-mass of the system. If the state of the sys- 
tem is well approximated by Eq. (23), the observed location 
would vary shot to shot with a probability distribution given 
by \4'a-n{xc)?' 0^ exp(-A^')'x^) (combined with any experimen- 
tal uncertainties associated with density measurement). For 



y < Jn 



this distribution would be wider than one would 



expect using the product approximation, most notably if the 
center-of-mass wavefunction is wider than a classical soliton 
width 1/ ^|yN ^ 1, this jumping effect would be most clearly 
visible. Non zero temperature would further increase this ef- 
fect, by introducing a statistical mixture of excited states of 
the center of mass. As a purely mechanical analogy, one could 
think of taking a photo of a swinging pendulum at a random 
time, the shape always looks the same but its position appears 
random. 
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VI. CONCLUSIONS 
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We study a ID system of identical Bosons with attrac- 
tive contact interactions, a Lieb-Liniger(-McGuire) gas, in the 
presence of a harmonic trapping potential. We present vari- 
ational and numerical many-body calculations, in both cases 
making use of the separability of the center-of-mass Hamilto- 
nian to split the problem into relative and center-of-mass de- 
grees of freedom. We use a unit system such that the Hartree 
soliton length is set to unity (ft = m - g{N - 1) = 1), leav- 
ing two parameters, the number of atoms, A^, and -\fy, the 
dimensionless ratio between the Hartree soliton and harmonic 
oscillator lengths. 

Our key results are firstly that we have derived a first order 
energy correction to the ground state of the relative degrees 
of freedom from the introduction of a harmonic oscillator po- 
tential [given in Eq. (27)], which is used in a variational mini- 
mization technique. This is proportional to y^ and the correc- 
tion per atom tends to the mean field prediction from below, 
the relative difference is less than 1% for A^ > 165. 

Secondly we have determined the validity range of y of our 
many body ansatz, consisting of the free many-body ground 
state with a Gaussian envelope as given in Eq. (23). Essen- 
tially as the trapped ground state deviates from this it becomes 
less "soliton like", we quantify this with the "effective single 
particle overlap", given by the A^th root of the overlap between 
a variationally obtained ground state and our ansatz. For A^ 
large, this overlap is greater than 0.99 for y < 0.16. Numerical 
calculations of energy in the strongly trapped region, y > 1, 
indicate energies are still considerably lower than the non in- 
teracting case. 

Thirdly we show, via a numerical investigation of over- 
lap between the free Hartree product solution and the free 
many-body ground state with a Gaussian envelope [given in 
Eq. (23)] describing the center-of-mass wavefunction, that the 
two wavefunctions can have high agreement, even at a many- 
body level. This high overlap occurs when the modulus square 
of the envelope function matches the probability distribution, 
associated with the Hartree product, for the center-of-mass 
position, which occurs when y » 0.16. However, current ex- 
periments with bright matter-wave solitons are such that the 
center of mass is localized to much less than a soliton width, 
indicating this is unlikely to be an observable effect. 

In addition to these physical results, we outlined a numeri- 
cal method for computing many body eigenstates, this uses a 
basis set of harmonic oscillator eigenstates, truncated at a par- 
ticular energy. This is then project a into a subspace of states 
with the center-of-mass wavefunction in a specific state, using 
the ladder operator for center-of-mass excitation. This makes 
use of the separability and achieves a reduction in the size of 
the basis set required by a factor of n/ V6-Ecut (where E^ut is a 
cut-off" energy) or better, greatly improving speed of the diag- 
onalization and allowing us to investigate the internal degrees 
of freedom separately. 
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discussions, and the UK EPSRC for funding (Grant No. 
EP/G056781/1) and the Jack Dodd Centre (S.A.G.) for sup- 
port. 



Appendix A: The ID free system and the full eigenspectrum via 
the Bethe Ansatz 



We briefly recapitulate aspects of the treatment of the at- 
tractively interacting Lieb-Liniger gas [ , ] in order to set 
notation within our chosen system of units. 

In order to find the solution to the ground state of Eq. (6) 
with y = 0, we note that the wavefunction in the region xj < 
X2 < ■ ■ ■ < xn is solved by 



H^ 



1 



Va^ 



z 

IP) 



A(P)sxp\iJ] 



, P'P{k)Xk 



k=\ 



(Al) 



where a sum over all permutations of the set !P = { 1 , . . . , A^) 
is performed to make it symmetric, the energy eigenvalue is 
thus simply equal to £ - YjkPlf^- Each permutation has a 
coefficient associated with it that is linked to the boundary 
conditions when x^ = xj.+i, for an interacting system they can 
be determined by the equation [37] 



A(P') = A(P) 



PP(k+i)- Pp(k)-i/(N - 1) 

pp(k+i) - Pnk) + '/(A^ - 1) ' 



(A2) 



where T" is the permutation swapping the k\h and {k + l)th in- 
dices, the coefficient of the identity permutation is determined 
by the normalization condition. The center-of-mass motion is 
independent of the interactions, and so will have eigenstates 
of plane waves. In the case of attractive interactions. These 
momenta can also have very specific imaginary components 
corresponding to bound clusters of atoms. The ground state 
of relative motion occurs for 



Pk 



.N +l-2k 
' 2{N - 1) 



(A3) 



in which all the permutation coefficients apart from one (the 
identity permutation) are equal to zero. Higher eigenstates 
can have multiple bound state clusters or strings, each with 
an associated real momentum P„, and imaginary components 
(that must sum to zero) which are spaced in units of 1 l{N - 
1). If we have 77 clusters, each of size n,„, we have momenta 
associated with the mth cluster given by 



Pk = Pm + i- 



2™ I n( + \-2k 



m-\ 



2{N-\) 






This 



state would 



^Pn)^ 



normally 
the total 



be 
energy 



ne. (A4) 
c=\ 

denoted as 
of the state. 



|«i,pi,n2,P2,-- 

E - YikP\l'^7 scales as though these are 77 isolated single 

soliton states and thus the energy eigenvalue is 



y hmPj 

117= 1 ' 



»/»(»,; - 1) 

24(A^ - 1)2 



(A5) 
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The total number of different combinations of clusters scales 
p{N), the number of ways to partition A^ with integers (cf. 
Appendix F4 and [48]) 



Appendix B: Normalization, energy and overlap using the 
variational state 



1. Preamble 



wave-function 



iM)\ 



exp 



( N k-\ \ 

V k=2 j=\ 

Ny 



ZXk_ 
N 






X exp 



dp 



exp{-p^/2y) 



7.nj 



N k-1 



ip 



ZXk _ " V^ V^ 
T7X7 TZjZj 



Va^ 



^k ^ i 



k=2 j=l 



(Bl) 



In order to make use of our variational state given in 
Eq. (31), we must calculate the normalization constant and 
expectation value of energy. Calculations for the energy and 
normalization constants for all the eigenstates in free space 
(y = 0) can readily be found in literature [ ], it is also the 
case that the center-of-mass component of the Hamiltonian 
can be considered separately and so taking a finite center-of- 
mass component does not significantly alter the calculations. 
However, the choice of normalization condition for a non local 
system is somewhat arbitrary and conventions vary between 
papers, we choose a normalization that means both the relative 
and center-of-mass parts are normalized to unity with respect 
to Jacobi coordinates. Also most derivations of the energy rely 
on the fact that the gradient discontinuity at the points Xf, - Xj 
in the wave-function, exactly cancel the interaction terms (es- 
sentially from the condition of being an eigenstate), and thus 
these terms can simply be ignored. Because of our variation 
of A, this will no longer be the case and thus we are forced 
to make a more explicit calculation of the kinetic energy. In 
addition to this we derive a first order energy correction to the 
relative degrees of freedom, which is a new result. 



with cr - 1/(N - 1) corresponding to Eq. (11), however 
for greater generality we allow this parameter to be free in 
order to use these results for variational calculations where 
cr — > AjiN - 1), which would correspond to Eq. (31). Also 
one may wish to consider instead units in which the harmonic 
oscillator frequency and length are set to unity and the in- 
teraction constant rescaled to g, in which case the replace- 
ment cr — > A\g\ would be used instead, or indeed in S.l. units 
cr — > AM\gid\/}f. Essentially this term serves to allow easy 
conversion between unit systems and making variational ma- 
nipulation easier. 

In the form of Eq. (Bl), it is far simpler to perform the 
integrals of the coordinate variables. Calculating (i/'tarl'/'var) 
in coordinate space will require integration over A^ spatial in- 
tegrals and two momentum integrals, however we only need 
to integrate over the simplex region xi < X2 . . . < xai as by 
Bose symmetry any integration over any such region will be 
identical, hence we multiply by factor of A^! to include all pos- 
sibilities for such a regions construction. Within this simplex 
region, all arguments in the absolute value signs are positive 
and the wave-function is given by 



(xi,..,JCAf|lA™r) =M 






dp 



^wi-p /^y) 



Iny 






^/N 



(B2) 



2. Normalization 



In order to normalize Eq. (11), we could insist that two 
states with different center-of-mass momenta are orthonormal, 
i.e. {p',N\p,N) = 6{p' - p) such as was calculated in [37] 
or consider wave-function to be trapped in a box which we 
allow to grow infinitely large [ ]. However we are inter- 
ested in the normalization of the free space solution with a 
Gaussian center-of-mass envelope and freedom to tune a vari- 
ational parameter, denoted i^vai in Eq. (31). This result and 
technique will also be used in appendix B 4 and follows the 
method of [ ] . We consider a Fourier decomposition of the 



with I3{k) - (N + I - 2k)cr, and using the notation 
f = r^ dxi f'' dxo . . . r dxfj, we can now 

J—CC<X^<X-)<...<X^<CO J— OO J—CC J—CC 

express the inner product as 



*t:vs>=«wi /£>.,<<.= ?^'^-"''*^'>''^' 



X 



03<X\<...<X^<00 



Iny 
exp 2j — + fi(k)xk 



y/N 



(B3) 



transformation of variables p - ipi+ p2)/'2 and p' - Pi— pi 
with Jacobian unity then allows us to perform the integral over 
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p leaving 



(-A^lll-AllJ-) = N\N^ 



dp' 

oo 



exp(-p'^/4y) 



X 



-cxd<j:i<---<a,v<oo 



2^JJFy 



Xk 



(B4) 



To perform the remaining integrals we note that 
J_^ dxexpiax + by) = exp[(a + b)y]/a, denoting 



a(k) ^ Yj/3(1) ^ crk(k - N) , 



(B5) 



1=1 



and noting fl(A^) = 0, we can recursively use the previous re- 
sult to perform all but one of the spatial integrals and obtain 



dp 

oo 



, exp(-p'V4y)A(JV,pO 
2^lny 

I dxj^ exp iip'xj^ yfNj , (B6) 



with 



{-I 



An;,p') = Y] 



k=l 



a{k) + 



ip'k 



(B7) 



Integrating the final term gives 2nS(p' y/N) and the momen- 
tum integral is then trivial, noting that A(N,Q) - l/iN - 
l)!^cr^'^ '^ gives us the final result for the normalization factor 



N, 




l)\o-(N-i) 



(B8) 



As was mentioned before, the normalization factor for 
Eq. (31) in our units is obtained by letting cr = /t/(A^ - 1), 
in the case of A = I where this relates the ground state in in- 
finitesimal trapping Eq. (11), we refer to this constant simply 
as N. It is also worth noting that both the com wave-function 
and the relative are both chosen to be normalized to unity with 
respect to Jacobi coordinates, hence the (y/A^Tr)'^"* relates to 
the center-of-mass part and the rest to the relative component. 



3. Kinetic and Interaction energy 

We wish calculate the kinetic energy and potential energy of 
the variational state, i.e. the expectation of Eq. (6) with A set to 
zero on Eq. (31), which we will denote Hfiee. Due to the sep- 
arability of the wave-function and Hamiltonian, it is sufficient 
to consider only the relative part of the wave-function and note 
that the center-of-mass kinetic energy is given by j/A. We 
first denote ^(xi, ..,XAr) = expf-cr^^jZiti l-f* - '*^;l/2), be- 
ing the relative part of the variational wave-function (up to a 
normalization factor) and calculate the second derivative with 



respect to some coordinate X{ 



1^2 

2 5jc' 

cr 

'2 



^(f{Xu..,Xff) 



8 X 8^ 



dxe 



8x^ 



X(p(xi,..,xm) 
( 
^sgn(xf-jci) +2 ^6{x(-Xk) 



\kU 

X(p(xi,..,xn) 



( 



(B9) 



The first term in Eq. (B9) can be split up into terms of the 
form sgn2(xf - x/,) - 1, of which there are (A^ - 1) and terms 
of the form sgn(xf - Xa)sgn(x( - Xb) with a ^ /?, of which 
there are (A^ - 1)(A^ - 2). The former will evaluate to unity by 
normalization of the wave-function, however the latter terms 
will equal + 1 when x/ < x^ < x/, or X( < xi, < x^ and when 
Xa < Xh < X[ or Xh < Xa < X( and -1 when x^ < Xf < x^ or 
Xb < xt < xi,; the wave-function must be identical in all these 
6 simplicies due to Rose symmetry and so the expected value 
of these terms will equal 1/3. When the sum over all { is per- 
formed, these terms will total to A^(A^- l)cr2(l 4-(A^-2)/3)/4. 
The latter terms of the form 6{xc - Xk) can then be combined 
with those from the interaction part of the Hamiltonian, not- 
ing that there are twice as many terms but 6{a - b) - 6{b - a). 
Reinstating cr - A/{N - 1) we have 



N k-l 



<<At'il 



k=l k k=2 1=1 



A^N 



8(A?-1)\ "^ 3 ) "^4 



k=2 ,/=l 

r 



^-1 

N- 1 



N k-l 



k=2 7=1 



(BIO) 



all that remains now to calculate the value of the expecta- 
tion value of the delta function terms. Following the method 
in Appendix B 2 we integrate over a simplex region -oo < 
X[ < X2... < xn < oo, as a result of this we need only 
consider the N - 1 terms of the form S(xk - Xk+i) as the 
rest will be zero. Each integral will be the same as in Ap- 
pendix B 2 except missing a factor of 2/a{k) for each term 
6(xk - Xk+i), hence the result will equal XfJi' a(k)/2 (using 
the result J ^ dxf(x,y)6(x - y) - f(y,y)/2). Hence 



N k-l N-1 

k=2 i=\ k=\ 

_ A{N +\)N{N - 1) 

" \2{N-\) 



, (BID 
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finally, substituting in this result into Eq. (BIO) we have 



N k-l 



<^vi;-i4z|?-]^ZZ'^(---M^> 



t=i ^""k 



k=2 j=l 



N(N + 1) , r 

— -[-A^ + 2A(A - 1)] + - 

24(N - 1) 4 

Eo(2A -A^) + ^, 



(B12) 



where £0 = -A^(A^ + l)/24(N - 1). This discussion has 
not mentioned the harmonic envelope of the center-of-mass 
function, however due to the separability of the Hamiltonian 
this will only add factor of the center-of-mass kinetic energy, 
which will be independent of A regardless of what the center- 
of-mass wave-function is. This energy term is combined with 
the potential energy calculation derived in Appendix B 4 to 
form the basis for a variational principle. 



4. Derivation of the first-order perturbation energy 



,(i)i 



,M)\ 



This section is related to the calculation of (lAvarl^WliAvar). 
where V(x) - y^ Ylk^\ -^^/^^ note that Eq. (27) is given by this 
quantity minus the center-of-mass energy. It is again easier 
not to perform this integral in Jacobi coordinates but to Fourier 
transform out the center-of-mass, we will also again replace 
the factor 1/(A^ - 1) with a to generalize the results for our 
variational principle. This calculation is similar to, although 
more complicated than, the calculation of the normalization 
factor; to that end we can start the calculation from Eq. (B4) 
(as no spatial integrals are yet performed) adding in the poten- 
tial factor Y{x) giving 



2tiA{N, 0) J_o„ 



4y 



X 



oo<.Vl <.V2<...<Xw<oo 



Z- 



X exp ^ ip — ^ + y6(fc)xi 



Va^ 



(B13) 



The same recursive integral procedure can be applied here ex- 
cept we need two additional results, true for real(^) > 



xj —Oi. 



{\ 



dx x" exp(A:x) 



exp(^y) 
^ exp(^y) 

(kyf-lky+l 



ifn = 0, 

if« = l. (B14) 



k^ 



exp(A:y) if n = 2. 



Let us consider only the latter part of Eq. (B 1 3) omitting the 
constant ^jNy^ IAnA{N,0), taking the integrals in order from 
xi to xn, the integral over X( will be over a function of the 
form 



m^A({,p')\kQ(() + h{€)xt + k2(()x]+ 2 x], 



e'=c+\ 



X exp a({)xe, + ^ /3i{')xe' 



('={+1 



ip{xc 



(B15) 



with kail) = ki{l) = and kzil) = 1 and A{(,p') defined in 
Eq. (B7). The common prefactor of A{{,p') is the equivalent 
of k from Eq. (B14). Besides this, each integral will increase 
the factor in front of the x^ term by one each time and hence 
kiH) - I. Contributions to A;i(^ -H 1) come from ki{€) and k2{() 
and as such Eq. (B14) impUes ki{{+ 1) = ki{()-2kz{(){a{i) + 
ipti yfNY^, given that ki{\) - this impHes 



k,(,(+l)^-2Y^ 



^-f a{k) + 



ikp' 



^i(A^)lp'=o 



N-l 



(B16) 
(B17) 



Applying this same induction logic down to ko gives 

ki({) ^ k2(i) 



ko({ + 1) = koit) - 



a(€) + 



ip't 



+ 2 



(«W+^)' 



f e 
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\fN' 
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N-l t 
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(B18) 



simply evaluating these at A^ and performing the final integra- 
tion over xm then yields 



.7 -a 



dxt^ I{N) - 



2A(N, p')7T 
yfN 



k2(N)6"(p') 
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(B19) 



we then insert this expression back into Eq. (B 13) giving 



<iA<li|V(x)|(^S.) 



r 



Lj^'^H^ 



2A{N,Q)J_^ ■ "\ 4r 
k2(N)6"(p') .ki(N)6'(p') 



N 



+ i- 



^fN 



■ ko(N)6(p') 
(B20) 



The integral over the 6(p') term can be performed immedi- 
ately and gives y^kQ(N)l2. Considering next the integral over 
6'{p'y, since exp(-;?'^/47') has zero gradient at the origin it 
will not contribute, however the terms 



d 2/ 

^^i(A^)l;y=o=— ^^,2 



N-l ^ 



-A{N, 



N-l J 

,/)lp'=() = ^A(A^,0)y-, 

o-^/N f^k 



dp o-y/N 

will contribute to Eq. (B19), giving 



(B21) 
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(B22) 
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Finally for the 6"(p') term we must include A"{N,p')\p'=o - 
-AiN,0) i;^fii(l + 5ki)k{INa(k)a(() and the differential of a 
Gaussian, hence we have 






dp'6"{p')e\p 



-P 
4y 



A(l,p') 



■■A(N,0) 



1 1 

2^ "^ N^ 



N-\ , Y N-l - 



4=1 



*=1 



k^ 



(B23) 



Summing these three terms together, and substituting k2(N) 
N, we are left with 



5. Overlap of the relative components of the variational 
wavefunctions 

Finally we consider the overlap between the relative parts 
of the variational wavefunction with A > I and the ground 
state in infinitesimal trapping 7 = 0, A = I, given by 



((Atlil-Aw) =M/^^ I dxi 



dXN\i/a 



X exp 



dx\ ■ ■ ■ I 

oo «7 — ex 



^ ^+1 " '-' 



2{N-\) 
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(B29) 



(•A^^IVI^tli) =^ 
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The first term in this expression is equal to 7/4, which is sim- 
ply the potential energy of the center-of-mass component.lt 
can be proved via induction [49] that the double sum is equal 
to 



N-l 1^1 1 IfN-l , Y N-l T 



,^_,m-i){^,N-k 2N[f^^k^ 



k=l 



, (B25) 



thus reinstating cr - A/(N- 1), the remaining terms simplified 
down to 



<^t'ilV|^l,li) = 






yHN-lf 



NA^ 



k^l 



k^ 4 ' 



(B26) 



which is used in Sec. 111. Equation Eq. (27) is the first order 
energy correction to the free soliton with Gaussian center-of- 
mass envelope and is obtained by subtracting the center-of- 
mass energy and setting A - I 



£« = 



j^(N - l) 
N~ 



2 N-l 



§^^' 



(B27) 



a. Energy correction from potentials of higher powers ofx 

An energy correction for general power law potentials can 
be derived in the mean field case. For Re(m) > -1 



f 



dx 



sechix/iy \xr _ I m\((m)(l - 2'-'") m ?t 1 , 
4 2~"1log(2) 



m — I , 
(B28) 
with ({m) the Riemann zeta function. A similar result would 
be desirable to calculate energy correction from an anhar- 
monic potential for a quantum soliton, although potentials 
with m ^ 2, will couple the center-of-mass and relative de- 
grees of freedom together (possibly only very weakly) and so 
[Hem, Hys\] + and this is only of limited use. 



This calculation can be achieved by performing the calcula- 
tions in Appendix B 2 with cr — > (1 -h A)I2{N - 1), the re- 
sulting factor will not equal unity and instead will be equal to 
A/^i Ar^/Ar?j^ (wj. Therefore that the overlap is given by 



(-AU'rl-A^i) = 



^(A'-l)/2 



ntY(i + ^)/2 



2^[A 
(l+^)j 



(B30) 



which is used in Sec 111 C. 



Appendix C: Ground state energy for //r(^2) in the interaction 
dominated regime 



Using the identity [50] 

r(z+i/2) 



r(z) 



k=2 ^ ) 



(CI) 



where the c^ are coefficients for the higher order terms in the 
asymptotic expansion, we see from Eq. (20) that 



/-V{) 



1 H- 



1 °° 



Ck 



{-y^f 



1 



2V2^ 



(C2) 



Hence, taking the limit y — > (interaction dominated 
regime) implies vq — > -oo , and we may truncate the asymp- 
totic series. To lowest order y-vo * 1/2 -\/2y, which we sub- 
stitute into the right hand side of [rearranged from Eq. (C2)] 



-vo = ^ + ;^+0(Von 



272^ 8V=^ 
squaring the result to get 



(C3) 



vo = -- 



\-\-oe) 



(C4) 



Hence, substituting Eq. (C4) into Eq. (21) for n = yields 



lim£R,o = --+0(/). 

•)'->0 4 



(C5) 



17 



Appendix D: Energy correction to the Hartree product state 

This section derives the energy correction to the Hartree 
product state \^h}, this state is a product of A^ identical single 
particle wavefunction <l>(jc) = sech(ji:/2)/2. The first result 
we require is the potential energy correction to each single 
particle wavefunction 



1 r°° 

■^ ^— CO 



r 



r 



dx x^sech^(x/2) 



2 2 



(Dl) 



the total energy correction is thus A^ times this value. However 
we are interested only in the relative energy correction given 
by 






r (-AhI 



.A'-l 
IN 






XkXj\iJ/}i} ■ (D2) 



Appendix F: Truncating the Hilbert space by introducing 

energy cut-offs and projection to the zero center-of-mass 

excitation subspace 

1. Integer partition and energy level degeneracy 

The relation between energy level degeneracy in systems 
of identical particles and number partitioning has been inves- 
tigated in [51-53] and references therein. Within a one di- 
mensional Harmonic oscillator, the degeneracy for A^ distin- 
guishable particles scales the same as the degeneracy for one 
particle in an A^ dimensional spherically symmetric potential. 
This is not the case for indistinguishable particles, to calculate 
these we must use introduce integer partition functions. We 
introduce the notation p([a,b],m) being the number of ways 
to partition an integer m using only integers a < z < b, in order 
to compute these for a given b, we use the recursion relation 



p([a,b],m) 



k=\ 



if a > min(m, b) and m + 

1 if [fl = m or m-Q\8i.a<b 

p([a + l,b],m) + p([a, b],m - a) otherwise . 

(Fl) 



All the cross terms of the form x^Xj will evaluate to zero as 
sech(x) is an even function, thus leaving only the power terms. 
By Bose symmetry (/(x;t)) = {/(xj)} and thus the value of all 
the terms in the first sum will be identical to the single particle 
correction and we have 



4" = (A^-i)^ 



(D3) 



Appendix E: Exact solution to the variational minimization 

The solution derived to the minimization equation (35) is 
given by 



A\A -l)-K^O, 



(El) 



with *■ > defined by Eq. (36). This equation has exactly one 
real positive solution Aq corresponding to an energy minimum, 
this solution can be derived analytically [43] (cf. Ref. [30]); it 
is given by: 



io = - (l -H Va -H V3-Ah-2A-'/2) 



with 



A= 1 



Y \3 



2/3 



Y = (-9;^ + V3 V27^?T256?) 
A Taylor expansion about a: = yields 



1/3 



Ao^I+k-3k^ + 0{k^). 



(E2) 

(E3) 
(E4) 

(E5) 



This works by noting that we can divide a partition into two 
distinct sets, partitions which uses only numbers larger than 
fl, being p{{a + \,b],m), and partitions which uses a at least 
once in the partitions, p{{a, b],m- a). Also p([a, ^], 0) = by 
convention. 

Using the usual Fock space representation of these har- 
monic oscillator states |A^o, N], . . .} with Yjk^k = N and defin- 
ing E as the energy of the state (with no interactions) minus 
the ground state energy divided by y 



E 

7 



N 
2 



Y^kNk 



(F2) 



k=0 



Given that each occupancy of the A:th mode raises the energy 
by k it can be seen that the degeneracy of the energy level E is 
given by the number of ways to partition E using A^ nonnega- 
tive integers. Denoting (t>(E, €) as the ways to partition E ini 
numbers we have 



g{E,N) ^Yj^{EJ) 

(=0 



(F3) 



It is also known that this sum is equal to the number of ways to 
partition an integer '£' using only numbers less than or equal 
to N i.e. g(E,N) = p([l,N],E). 



2. Truncation with an energy cut-off 

In order to make a basis computationally manageable, it 
must truncated to be made to be finite. This is achieved by 
only taking states with energy less than an arbitrary cut-ofF rj, 
note that this also implies that A^^ = if ^ > rj. The size of 
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this truncated Fock state basis is given by 



Y^giE,N) = J]pi[l,N],E). 



(F4) 



£=0 



£=0 



The reason an energy cut-ofF is chosen rather than a mode cut- 
off at T] (ahhough as mentioned before this is impHcit in an 
energy cut-off method) is two fold. Firstly in order to project 
into the center-of-mass and relative excitation basis we require 
all the states with a given energy E [the Hamiltonian (52) is 
block diagonal], if we do not have all those states the pro- 
jection is not possible. Secondly having just a mode cut-ofF 
would include the state |0, 0, . . . , A^) with E = Nrj, but not the 
state |A^ - 1, 0, ... 0, 1) (one occupancy in the ;; + 1th mode) 
with £■ = 77 H- 1, as long as harmonic oscillator energy remains 
a non negligible quantity, the former state will have almost no 
mixing to the ground state, making it a very inefficient trunca- 
tion. 



3. Deriving the projector to the center-of-mass basis 

As we have expressed the Hamiltonian (48) in terms of a^ 
and a^., the creation and annihilation operator for bosons in 
mode k, it is far simpler to compute the matrix elements in 
terms of basis states in the \No,Ni, . . .), occupation notation. 
Therefore we wish to calculate the elements and then project 
into eigenstates of the center-of-mass Hamiltonian H^m given 
in Eq. (52). It is therefore sufficient to diagonalize A^A^, 
[given by Eq. (5 1 )] as this is the only operator dependence in 
Hem, using basis states of the form \N{),Nu . . .). This gives a 
square matrix P of eigenvectors of center-of-mass, which can 
project the truncated Fock state basis into this new basis, and a 
vector of eigenvalues. This is computationally simple as A'^A'' 
cannot mix states of different energies and therefore is block 
diagonal when states are ordered by energy and each block can 
be diagonalized separately. By removing all the columns of P 
with associated eigenvalues not equal to zero (meaning they 
have excitations in the center-of-mass mode) we are left with 
a rectangular matrix P which projects into ffiis ground state of 
center-of-mass excitation subspace that we call the 'reduced 
basis'. 

Using P results in a far smaller basis set (discussed in the 
next subsection) without changing any of the relative dynam- 
ics, however it is not immediately clear what states in this 
new basis relate to. Given that each partition of E into A^ pos- 
itive integers has the interpretation that each integer k repre- 
sents a single occupancy in the A:th mode, one may ask what 
the relation to quantum numbers is of partitions in terms of 
integers less than or equal to A^, for instance E - 2 can be 
partitioned by 1 -1- 1 and 2. Given that we know the ladder op- 
erator associated with the center-of-mass mode A^ , satisfies 
[Hq,A'^~\ - +yA'^ and is thus spaced in steps of unity times y, 
we can associate all the 1 's in a given partition with a quanta 
in this mode. Assuming we have ( quanta in the center-of- 
mass mode, this leaves all the numbers 2 < z < A^ as ways to 
partition E - H, which must then relate to some relative exci- 
tation modes. Going back to £ = 2 the partition, 2 = 1-1-1 



is two quanta in the center-of-mass mode i.e. A'^A'^\N,Q, . . .) 
and the partition 2 = 2 is one quanta in the ffist relative mode. 
In order to help understand this we examine the N -2 case 
in ffist quantization, using Jacobi coordinates [Eq. (15)]. The 
Hamiltonian can be expressed in two commuting parts 



H, 



rel 



1 52 2 

'4 5x2^''^' 






and 



(F5) 
(F6) 



For distinguishable atoms, these would each have normal har- 
monic oscillator eigenstates (up to a scaling factor), which 
can be multiplied together to create a many-body eigenstate. 
However, we require Bose symmetry of the many-body wave- 
function: il/(xi,X2) = ij/(x2,xiy, in terms of Jacobi coordi- 
nates this implies no conditions on xq but that iff(xc,^2) = 
il)(xc,-^2) and hence odd eigenstates for 7/^1 are disallowed 
and relative energy levels are spaced in units of 2. 



4. Basis size reduction 

As mentioned in Appendix F3, the center-of-mass mode 
ladder operator A'^ of Eq. (51) has an energy spacing of 
unity, implying relative excitation modes are spaced in units 
of 2,3,.., A?. 

In our reduced basis, the subset with the center-of-mass in 
the ground state, we can no longer partition E using the num- 
ber 1 . Therefore the energy degeneracy g{E) of level E in 
the reduced basis, is ffie number of ways to partition E using 
integers z satisfying 2 < z < N, i.e. g(E,N) = p([2,N],E). 
Therefore the number of basis states in the reduced basis rel- 
ative to the occupation number basis with cut of E is given 
by 



A(/7,A^) = 



ZI_^P([2,N],E) 



■'£=0 



EtoM[l,A^],£) 



(F7) 



We can use the equation of Eq. (Fl) to write p{[2,N],E) = 
p{[l,N],E) - p{[\,N],E - 1), in ffie sum from to rj, all 
terms cancel apart from those at the end points of the sum, 
leaving only the term p([l,A^],77) and hence the size of the 
reduced basis is just the degeneracy of the //th energy level in 
the occupation number basis, thus we have 



A(t],N) 



PJlhNlrj) 
Zl^^p([l,N],E) 



(F8) 



Essentially this property can be seen from projecting the set 
of kets with energy rj into the center-of-mass excitation basis, 
this set will contain all the relative excited states with energy 
less than or equal to 77, but with additional center-of-mass ex- 
citation. 

The basis reduction for N - 2 can be calculated by noting 
there are [kl2\ + 1 ways to partition k using 1 and 2 (the nota- 
tion [k\ means round k down to an integer), thus the reduced 
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FIG. 5. (Color online): Reduced basis size divided by truncated basis 
size, given by Eq. (F8), for different cut-off energies. Top to bottom 
lines are for cut-off energies t; = 10, 20, 40, 60, 80, 100, dotted lines 
are the estimate of Eq. (Fl 1). Basis reduction is most significant for 
small N but Eq. (Fl 1) provides a good upper bound on reduction for 
large N. 



basis is [rj/l] -)- 1 in size, the number of states in the truncated 



occupation number basis is 



2(L^/2J + 1) = ] 



1+77-1- '7^/4 



if T] even 



ri + (if - l)/4 if rj odd. 



(F9) 



To leading order the reduction A(77, A^) goes as 2///. Such sim- 
ple analytic expressions are not known for general A^, however 
we have the following expression by Ramanujan [ .ii] 



p{[l,N>n],rf) 



1 



4?7V3 



exp TT 




as rj 



(FIO) 



this can be used to get an asymptotic estimate of the basis 
reduction by replacing the sum in Eq. (F8) with an integral, 
giving 






^ 






(Fll) 



which will be our best estimate for the reduction achieved for 
large A^, note that this improves slower than the cc l/rj reduc- 
tion for the N - 2 case. This asymptotic estimate is included 
in Fig. 5, along with the reduction for intermediate values of 
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